clear all
*cap log close
set more off
set seed 603

* load utility programs --------------------
qui do "${code}/utility/ebayes"


* figure showing distribution of officer bunching propensities ------------------------

* panel a: raw propensity --------------------
use "${data}/out/4-main", clear
collapse (mean) discount, by(officerid)

#delimit ;
hist discount, 
lcolor(dknavy) fcolor(navy%10) gap(30) bin(60) 
frac plotregion(lcolor(black) lwidth(medium))
graphregion(color(white)) ylab(,nogrid) xlab(,nogrid)
xtitle("Bunching Propensity (Raw)") 
ytitle("Fraction of Officers") ;
#delimit cr
graph export "${out}/apx_iv/propensity_raw.pdf", replace

* panel b: residualized ----------------------
use "${data}/out/4-main", clear

reghdfe discount, absorb(totfe) nocons resid
predict rdiscount, resid
collapse (mean) rdiscount, by(officerid)

#delimit ;
hist rdiscount if abs(rdiscount)<0.6, 
lcolor(dknavy) fcolor(navy%10) gap(30) bin(60) 
frac plotregion(lcolor(black) lwidth(medium))
graphregion(color(white)) ylab(,nogrid) xlab(,nogrid)
xtitle("Bunching Propensity (Residualized)") 
ytitle("Fraction of Officers") ;
#delimit cr
graph export "${out}/apx_iv/propensity_resid.pdf", replace

* panel c: empirical bayes ---------------------
use "${est}/getfe", clear
summ est, d
replace est = est-r(p50)
ebayes est se, gen(est_eb)

kdensity est,    gen(x0 y0) nodraw
kdensity est_eb, gen(x1 y1) nodraw

summ est
local sig1 = "Raw: {&sigma}{superscript:2} = `:di %4.3f `r(sd)'^2'"
summ est_eb
local sig2 = "EB: {&sigma}{superscript:2} = `:di %4.3f `r(sd)'^2'"

#delimit ;
line y0 x0, lwidth(medium) lcolor(dknavy) ||
line y1 x1, lwidth(medium) lcolor(dkgreen) lpattern(dash)
graphregion(color(white)) plotregion(lcolor(black) lwidth(medium))
ylab(0(0.5)2.5,nogrid) xlab(,nogrid) ytitle("Density") xtitle("Estimated Officer Effect")
legend(pos(2) ring(0) cols(1) region(lstyle(border))
lab(1 "Raw Effects") lab(2 "Empirical Bayes"))
text(1.5  0.98 "`sig1'", place(w) size(med)) 
text(1.36 0.98 "`sig2'", place(w) size(med)) ;
#delimit cr
graph export "${out}/apx_iv/propensity_eb.pdf", replace
* ---------------------------------------------------------------------------------------


* figure showing correlation in propensities --------------------------------------------
#delimit ;
global cov = "female age agesq age_miss race_b race_h race_o race_u priorprison local
logzip zipincome_miss logprice veh_miss speed_py1 other_py1 crashany_py1" ;
#delimit cr

* panel a: locations -----------------------
use "${data}/out/4-main", clear
gen bunch = (speedd==9)
gen partition = partCOUNTY 
gen N = 1 
egen jidpart  = group(officerid partition)

reghdfe bunch ${cov}, absorb(fe=jidpart totfe, savefe) resid

collapse (mean) fe (sum) N, by(officerid partition)
drop if mi(officerid) | mi(partition )
reshape wide fe N, i(officerid) j(partition)
keep if !mi(fe1)&!mi(fe2)
gen Nw = N1+N2

/// Correlation estimates -------------
reg fe2 fe1, r
local bprint1 = "Unweighted: {&beta} = `:di %4.3f _b[fe1]' (`:di %4.3f _se[fe1]')"
reg fe2 fe1 [aw=Nw], r
local bprint2 = "Weighted: {&beta} = `:di %4.3f _b[fe1]' (`:di %4.3f _se[fe1]')"

/// Build plot -------------------------
#delimit ;
scatter fe2 fe1 [aw=Nw], mcolor(emidblue) msymbol(Oh) msize(tiny) ||
lfit fe2 fe1 [aw=Nw], lcolor(forest_green) lwidth(medthick)  ||
lfit fe1 fe1 [aw=Nw], lcolor(maroon) lpattern(dash)
graphregion(color(white)) plotregion(lcolor(black) lwidth(medium)) 
ylab(-1(0.5)1,nogrid) xlab(-1(0.5)1,nogrid) legend(off) 
xtitle("First Partition") ytitle("Second Partition")
text(0.9 -0.95 "`bprint1'" "`bprint2'", place(se)) ;
#delimit cr
graph export "${out}/apx_iv/corrs_location.pdf", replace


* panel b: time periods -------------------
use "${data}/out/4-main", clear

// Build partitions -----------
gen date = offensedate
bysort officerid: egen middate = median(date)
gen partition = (date>=middate)+1

gen bunch = (speedd==9)
gen N = 1 
egen jidpart  = group(officerid partition)
reghdfe bunch ${cov}, absorb(fe=jidpart totfe, savefe) resid

collapse (mean) fe (sum) N, by(officerid partition)
drop if mi(officerid) | mi(partition )
reshape wide fe N, i(officerid) j(partition)
keep if !mi(fe1)&!mi(fe2)
gen Nw = N1+N2

/// Correlation estimates -------------
reg fe2 fe1, r
local bprint1 = "Unweighted: {&beta} = `:di %4.3f _b[fe1]' (`:di %4.3f _se[fe1]')"
reg fe2 fe1 [aw=Nw], r
local bprint2 = "Weighted: {&beta} = `:di %4.3f _b[fe1]' (`:di %4.3f _se[fe1]')"

/// Build plot -------------------------
#delimit ;
scatter fe2 fe1 [aw=Nw], mcolor(emidblue) msymbol(Oh) msize(tiny) ||
lfit fe2 fe1 [aw=Nw], lcolor(forest_green) lwidth(medthick)  ||
lfit fe1 fe1 [aw=Nw], lcolor(maroon) lpattern(dash)
graphregion(color(white)) plotregion(lcolor(black) lwidth(medium)) 
ylab(-1(0.5)1,nogrid) xlab(-1(0.5)1,nogrid) legend(off)
xtitle("First Partition") ytitle("Second Partition") 
text(0.9 -0.95 "`bprint1'" "`bprint2'", place(se)) ;
#delimit cr
graph export "${out}/apx_iv/corrs_time.pdf", replace
* ---------------------------------------------------------------------------------------


* experience profiles (mentioned in text) -----------------------------------------------
use "${data}/out/3-officer", clear
gen startdate = firstfhpft
replace startdate = firstfhp if mi(startdate)
keep officerid startdate
tempfile startdate
save    `startdate'

use "${data}/out/4-main", clear
merge m:1 officerid using `startdate', keep(3) nogen
gen expdays = offensedate - startdate 
summ expdays, d
replace expdays = r(p99) if expdays > r(p99)

gen exp = int(expdays/(365/12))
label var exp "Experience in Months"
gen expsq = exp^2

reghdfe discount exp expsq, absorb(totfe officerid) vce(cluster officerid)
test exp expsq
* ---------------------------------------------------------------------------------------


* relationship between propensity and officer characteristics ---------------------------
use "${data}/out/4-main", clear
gen rlenience = discount
reghdfe discount, absorb(totfe) nocons resid
predict alenience, resid
gen     alenience_sd = alenience
collapse (mean) rlenience (mean) alenience (sd) alenience_sd, by(officerid)
tempfile temp 
save    `temp'

use "${data}/out/4-main", clear
keep officerid lenient
duplicates drop
merge 1:1 officerid using `temp', keep(1 3) nogen
merge 1:1 officerid using "${data}/out/3-officer", keep(1 3) nogen
gen age  = 2007-yob
gen ageorig = age
replace age = age/10
gen age2 = age^2
label var age "Age"
label var age2 "Age Squared"
gen stemp = "01/01/2007"
gen temp  = date(stemp,"MDY")
drop exp*
gen exp  = int((temp - firstfhp)/365)
replace exp = 0 if exp < 0
gen exporig = exp
replace exp = exp/10
gen exp2 = exp^2 
label var exp "Experience"
label var exp2 "Experience Squared"

gen black = (race == "B")
label var black "Race $=$ Black"
gen hisp  = (race == "H")
label var hisp "Race $=$ Hispanic"
gen other = (race != "W" & race !="B")
label var other "Race $=$ Other"
label var college "Any College"

global CHAR  = "female black hisp other age age2 exp exp2 college"
global TABLE = "female black hisp other age exp college"


* write results as data ----------------------

* setup variables ------
forval i = 1/4 {
	gen b_`i' = .
	gen s_`i' = .
	gen obscount_`i' = .
}
gen obscount_mean = .
replace obscount_mean =_N

gen lab = ""
local k = 1
foreach y in $TABLE {
	local lab: variable label `y'
	replace lab = "`lab'" if _n == `k'
	local ++k	
}

* means of indep vars -------------------
gen xm = .
local k=1
foreach y in female black hisp other ageorig exporig college {
	summ `y'
	replace xm = r(mean) if _n == `k'
	local ++k
}

* binary lenience --------------------
local k = 1
foreach y in $TABLE {
	reg lenient $CHAR, robust 
	replace b_1  =  _b[`y']  if _n == `k'
	replace s_1  =  _se[`y'] if _n == `k'
	local ++k
}

qui reg lenient $CHAR, robust 
replace obscount_1 = e(N)

gen ym_1 = . 
qui summ lenient 
replace ym_1 = r(mean)
format %5.3f ym_1

* continuous, raw --------------------
local k = 1
foreach y in $TABLE {
	reg rlenience $CHAR, robust 
	replace b_2  =  _b[`y']  if _n == `k'
	replace s_2  =  _se[`y'] if _n == `k'
	local ++k
}

qui reg rlenience $CHAR, robust 
replace obscount_2 = e(N)

gen ym_2 = . 
qui summ rlenience 
replace ym_2 = r(mean)
format %5.3f ym_2

* continuous, adjusted --------------
local k = 1
foreach y in $TABLE {
	reg alenience $CHAR, robust 
	replace b_3  =  _b[`y']  if _n == `k'
	replace s_3  =  _se[`y'] if _n == `k'
	local ++k
}

qui reg alenience $CHAR, robust 
replace obscount_3 = e(N)

gen ym_3 = . 
qui summ alenience 
replace ym_3 = r(mean)
format %5.3f ym_3

* continuous, weighted --------------
gen wt = 1/((alenience_sd)^2)

local k = 1
foreach y in $TABLE {
	reg alenience $CHAR [aw=wt], robust 
	replace b_4  =  _b[`y']  if _n == `k'
	replace s_4  =  _se[`y'] if _n == `k'
	local ++k
}

qui reg alenience $CHAR [aw=wt], robust 
replace obscount_4 = e(N)

gen ym_4 = . 
qui summ rlenience [aw=wt]
replace ym_4 = r(mean)
format %5.3f ym_4

* format results as data ---------------
forval i=1/4 {
	gen writeb_`i'= ""
	gen writes_`i' = ""
	gen writeym_`i' = ""
}
gen writexm = ""
replace writexm = string(xm,  "%5.2f") if abs(xm) > 1
replace writexm = string(xm,  "%5.3f") if abs(xm) < 1 & abs(xm) > .09999
replace writexm = string(xm,  "%10.4f") if abs(xm) < .09999 & abs(xm) > .009999

forval i=1/4 {
	replace writeb_`i' = string(b_`i', "%5.3f") if abs(b_`i') < 100 & abs(b_`i') > .09999
	replace writes_`i' = "(" + string(s_`i', "%5.3f") + ")" if abs(s_`i') < 1 & abs(s_`i') > .09999
	replace writeym_`i' = string(ym_`i', "%5.3f") 
	replace writeb_`i' = string(b_`i', "%10.4f") if abs(b_`i')< .09999 & abs(b_`i') > .009999
	replace writes_`i' = "(" + string(s_`i', "%10.4f") + ")" if abs(s_`i') < .09999 & abs(s_`i') > .009999
	replace writeb_`i' = string(b_`i', "%10.5f") if abs(b_`i')< .009999 & abs(b_`i') > .0009999
	replace writes_`i' = "(" + string(s_`i', "%10.5f") + ")" if abs(s_`i') < .009999 & abs(s_`i') > .0009999
	replace writeb_`i' = string(b_`i', "%11.6f") if abs(b_`i')< .0009999 & abs(b_`i') > .00009999
	replace writes_`i' = "(" + string(s_`i', "%11.6f") + ")" if abs(s_`i') < .0009999 & abs(s_`i') > .00009999
}

* post to LaTeX table -------------------
capture erase "${out}/apx_iv/table_officers.tex"
file open fh using "${out}/apx_iv/table_officers.tex", write replace

file write fh ///
	"&\multicolumn{1}{c}{} &\multicolumn{1}{c}{Binary}&\multicolumn{3}{c}{Continuous}    \\\cmidrule(lr){3-3}\cmidrule(lr){4-6}" _n ///
	"&\multicolumn{1}{c}{(1)}&\multicolumn{1}{c}{(2)}&\multicolumn{1}{c}{(3)}&\multicolumn{1}{c}{(4)}&\multicolumn{1}{c}{(5)}\\" _n ///
	"&\multicolumn{1}{c}{Mean}&\multicolumn{1}{c}{Lenient}&\multicolumn{1}{c}{Raw}&\multicolumn{1}{c}{Adjusted}&\multicolumn{1}{c}{Weighted}\\" _n ///
	"\hline" _n

local k=1
foreach y in $TABLE {
	local label = lab[`k']
	local m1 = writexm[`k']
	local m2 = writeb_1[`k']
	local m3 = writeb_2[`k']
	local m4 = writeb_3[`k']
	local m5 = writeb_4[`k']
	file write fh "`label'  & `m1' & `m2' & `m3' & `m4' & `m5' \\" _n

	local m1 = writes_1[`k']
	local m2 = writes_2[`k']
	local m3 = writes_3[`k']
	local m4 = writes_4[`k']
	file write fh "&  & `m1' & `m2' &`m3' & `m4' \\" _n
	file write fh "[1em]" _n 

	local ++k
}

* table footer -----------
file write fh "\hline" _n

local m1 = writeym_1 
local m2 = writeym_2 
local m3 = writeym_3 
local m4 = writeym_4 
file write fh "Mean  & --- & `m1' & `m2' & `m3' & `m4' \\" _n

local mm = obscount_mean
local m1 = obscount_1 
local m2 = obscount_2 
local m3 = obscount_3  
local m4 = obscount_4 
file write fh "Officers  & `mm' & `m1' & `m2' & `m3' & `m4' \\" _n

file close fh
macro drop fh
* ---------------------------------------------------------------------------------------


* show common support of stringency instrument ------------------------------------------
use "${data}/out/4-main", clear
reghdfe harsh Z, absorb(totfe, savefe) res
predict double pscore, xbd 
replace pscore = 0 if pscore<0 & !mi(pscore)
replace pscore = 1 if pscore>1 & !mi(pscore)


// Store info for printing ----
summ harsh if pscore==0
local p0 = "{it:Pr(D | p=0)} = `:di %4.3f `r(mean)''"
summ harsh if pscore==1 
local p1 = "{it:Pr(D | p=1)} = `:di %4.3f `r(mean)''"

#delimit ;
hist pscore if harsh==1, frac lcolor(dknavy) fcolor(none) bin(50) gap(20)
addplot(hist pscore if harsh==0, frac lcolor(none) fcolor(maroon%50) bin(50) gap(20) ||
pcarrowi 0.21 0.92 0.21 0.98 (9) "`p1'", color(black) ||
pcarrowi 0.21 0.08 0.21 0.02 (3) "`p0'" , color(black)  )
graphregion(color(white)) plotregion(lcolor(black) lwidth(medthin))
ylab(,nogrid) xlab(,nogrid)
xtitle("Propensity Score") ytitle("Fraction of Sample")
legend(ring(0) pos(11) order(1 2) region(lstyle(border))
lab(1 "Treated") lab(2 "Untreated"))  ;
#delimit cr
graph export "${out}/apx_extrap/support.pdf", replace
* ---------------------------------------------------------------------------------------


